﻿/*
* Copyright (c) 2010-2012 Tesla Engine Group
* 
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
* 
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
* 
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
* THE SOFTWARE.
*/


using System;
using System.Globalization;
using System.Runtime.InteropServices;
using Tesla.Util;

namespace Tesla.Math {

    /// <summary>
    /// Defines a 4x4 row-major Matrix (following DirectX's matrix layout). 
    /// Right handedness conventions used by default.
    /// </summary>
    [Serializable]
    [StructLayout(LayoutKind.Sequential)]
    public struct Matrix : IEquatable<Matrix> {

        /// <summary>
        /// Value at row 1, column 1 of the matrix
        /// </summary>
        public float M11;

        /// <summary>
        /// Value at row 1, column 2 of the matrix
        /// </summary>
        public float M12;

        /// <summary>
        /// Value at row 1, column 3 of the matrix
        /// </summary>
        public float M13;

        /// <summary>
        /// Value at row 1, column 4 of the matrix
        /// </summary>
        public float M14;

        /// <summary>
        /// Value at row 2, column 1 of the matrix
        /// </summary>
        public float M21;

        /// <summary>
        /// Value at row 2, column 2 of the matrix
        /// </summary>
        public float M22;

        /// <summary>
        /// Value at row 2, column 3 of the matrix
        /// </summary>
        public float M23;

        /// <summary>
        /// Value at row 2, column 4 of the matrix
        /// </summary>
        public float M24;

        /// <summary>
        /// Value at row 3, column 1 of the matrix
        /// </summary>
        public float M31;

        /// <summary>
        /// Value at row 3, column 2 of the matrix
        /// </summary>
        public float M32;

        /// <summary>
        /// Value at row 3, column 3 of the matrix
        /// </summary>
        public float M33;

        /// <summary>
        /// Value at row 3, column 4 of the matrix
        /// </summary>
        public float M34;

        /// <summary>
        /// Value at row 4, column 1 of the matrix
        /// </summary>
        public float M41;

        /// <summary>
        /// Value at row 4, column 2 of the matrix
        /// </summary>
        public float M42;

        /// <summary>
        /// Value at row 4, column 3 of the matrix
        /// </summary>
        public float M43;

        /// <summary>
        /// Value at row 4, column 4 of the matrix
        /// </summary>
        public float M44;

        private static Matrix _identity = new Matrix(1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 
            0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f);

        private static int _sizeInBytes = MemoryHelper.SizeOf<Matrix>();

        /// <summary>
        /// Get the identity matrix.
        /// </summary>
        public static Matrix Identity {
            get {
                return _identity;
            }
        }

        /// <summary>
        /// Gets or sets the translation vector (M41, M42, M43) of the matrix.
        /// </summary>
        public Vector3 Translation {
            get {
                Vector3 v;
                v.X = this.M41;
                v.Y = this.M42;
                v.Z = this.M43;
                return v;
            }
            set {
                this.M41 = value.X;
                this.M42 = value.Y;
                this.M43 = value.Z;
            }
        }

        /// <summary>
        /// Gets or sets the up vector (M21, M22, M23) of the matrix.
        /// </summary>
        public Vector3 Up {
            get {
                Vector3 v;
                v.X = this.M21;
                v.Y = this.M22;
                v.Z = this.M23;
                return v;
            }
            set {
                this.M21 = value.X;
                this.M22 = value.Y;
                this.M23 = value.Z;
            }
        }

        /// <summary>
        /// Gets or sets the down vector of the matrix (negation of Up vector).
        /// Setting will negate component values.
        /// </summary>
        public Vector3 Down {
            get {
                Vector3 v;
                v.X = -this.M21;
                v.Y = -this.M22;
                v.Z = -this.M23;
                return v;
            }
            set {
                this.M21 = -value.X;
                this.M22 = -value.Y;
                this.M23 = -value.Z;
            }
        }

        /// <summary>
        /// Gets or sets the backward vector (M31, M32, M33) of the matrix.
        /// </summary>
        public Vector3 Backward {
            get {
                Vector3 v;
                v.X = this.M31;
                v.Y = this.M32;
                v.Z = this.M33;
                return v;
            }
            set {
                this.M31 = value.X;
                this.M32 = value.Y;
                this.M33 = value.Z;
            }
        }

        /// <summary>
        /// Gets or sets the forward vector of the matrix (negation of backward vector).
        /// </summary>
        public Vector3 Forward {
            get {
                Vector3 v;
                v.X = -this.M31;
                v.Y = -this.M32;
                v.Z = -this.M33;
                return v;
            }
            set {
                this.M31 = -value.X;
                this.M32 = -value.Y;
                this.M33 = -value.Z;
            }
        }

        /// <summary>
        /// Gets or sets the right vector (M11, M12, M13) of the matrix.
        /// </summary>
        public Vector3 Right {
            get {
                Vector3 v;
                v.X = this.M11;
                v.Y = this.M12;
                v.Z = this.M13;
                return v;
            }
            set {
                this.M11 = value.X;
                this.M12 = value.Y;
                this.M13 = value.Z;
            }
        }

        /// <summary>
        /// Gets or sets the left vector of the matrix (negation of left vector).
        /// </summary>
        public Vector3 Left {
            get {
                Vector3 v;
                v.X = -this.M11;
                v.Y = -this.M12;
                v.Z = -this.M13;
                return v;
            }
            set {
                this.M11 = -value.X;
                this.M12 = -value.Y;
                this.M13 = -value.Z;
            }
        }

        /// <summary>
        /// Gets or sets an individual component based on its index
        /// in the matrix. The indices follow the named components, so
        /// the first element is located at (1,1) and the last element
        /// at (4,4). Indices outside this range return a value of zero.
        /// </summary>
        /// <param name="i">One-based row index valid between [1,4]</param>
        /// <param name="j">One-based column index valid between [1,4]</param>
        /// <returns>The matrix value</returns>
        public float this[int i, int j] {
            get {
                switch(i) {
                    case 1:
                        switch(j) {
                            case 1:
                                return M11;
                            case 2:
                                return M12;
                            case 3:
                                return M13;
                            case 4:
                                return M14;
                            default:
                                return 0;
                        }
                    case 2:
                        switch(j) {
                            case 1:
                                return M21;
                            case 2:
                                return M22;
                            case 3:
                                return M23;
                            case 4:
                                return M24;
                            default:
                                return 0;
                        }
                    case 3:
                        switch(j) {
                            case 1:
                                return M31;
                            case 2:
                                return M32;
                            case 3:
                                return M33;
                            case 4:
                                return M34;
                            default:
                                return 0;
                        }
                    case 4:
                        switch(j) {
                            case 1:
                                return M41;
                            case 2:
                                return M42;
                            case 3:
                                return M43;
                            case 4:
                                return M44;
                            default:
                                return 0;
                        }
                    default:
                        return 0;
                }
            }
            set {
                switch(i) {
                    case 1:
                        switch(j) {
                            case 1:
                                M11 = value;
                                break;
                            case 2:
                                M12 = value;
                                break;
                            case 3:
                                M13 = value;
                                break;
                            case 4:
                                M14 = value;
                                break;
                        }
                        break;
                    case 2:
                        switch(j) {
                            case 1:
                                M21 = value;
                                break;
                            case 2:
                                M22 = value;
                                break;
                            case 3:
                                M23 = value;
                                break;
                            case 4:
                                M24 = value;
                                break;
                        }
                        break;
                    case 3:
                        switch(j) {
                            case 1:
                                M31 = value;
                                break;
                            case 2:
                                M32 = value;
                                break;
                            case 3:
                                M33 = value;
                                break;
                            case 4:
                                M34 = value;
                                break;
                        }
                        break;
                    case 4:
                        switch(j) {
                            case 1:
                                M41 = value;
                                break;
                            case 2:
                                M42 = value;
                                break;
                            case 3:
                                M43 = value;
                                break;
                            case 4:
                                M44 = value;
                                break;
                        }
                        break;
                }
            }
        }

        /// <summary>
        /// Get the size of the Matrix structure in bytes.
        /// </summary>
        public static int SizeInBytes {
            get {
                return _sizeInBytes;
            }
        }

        /// <summary>
        /// Creates a new instance of Matrix
        /// </summary>
        /// <param name="m11">Element at row 1, column 1</param>
        /// <param name="m12">Element at row 1, column 2</param>
        /// <param name="m13">Element at row 1, column 3</param>
        /// <param name="m14">Element at row 1, column 4</param>
        /// <param name="m21">Element at row 2, column 1</param>
        /// <param name="m22">Element at row 2, column 2</param>
        /// <param name="m23">Element at row 2, column 3</param>
        /// <param name="m24">Element at row 2, column 4</param>
        /// <param name="m31">Element at row 3, column 1</param>
        /// <param name="m32">Element at row 3, column 2</param>
        /// <param name="m33">Element at row 3, column 3</param>
        /// <param name="m34">Element at row 3, column 4</param>
        /// <param name="m41">Element at row 4, column 1</param>
        /// <param name="m42">Element at row 4, column 2</param>
        /// <param name="m43">Element at row 4, column 3</param>
        /// <param name="m44">Element at row 4, column 4</param>
        public Matrix(float m11, float m12, float m13, float m14, float m21, float m22, float m23, float m24,
            float m31, float m32, float m33, float m34, float m41, float m42, float m43, float m44) {
                this.M11 = m11;
                this.M12 = m12;
                this.M13 = m13;
                this.M14 = m14;
                this.M21 = m21;
                this.M22 = m22;
                this.M23 = m23;
                this.M24 = m24;
                this.M31 = m31;
                this.M32 = m32;
                this.M33 = m33;
                this.M34 = m34;
                this.M41 = m41;
                this.M42 = m42;
                this.M43 = m43;
                this.M44 = m44;
        }

        /// <summary>
        /// Adds two matrices together.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>Added matrix</returns>
        public static Matrix Add(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = a.M11 + b.M11;
            m.M12 = a.M12 + b.M12;
            m.M13 = a.M13 + b.M13;
            m.M14 = a.M14 + b.M14;
            m.M21 = a.M21 + b.M21;
            m.M22 = a.M22 + b.M22;
            m.M23 = a.M23 + b.M23;
            m.M24 = a.M24 + b.M24;
            m.M31 = a.M31 + b.M31;
            m.M32 = a.M32 + b.M32;
            m.M33 = a.M33 + b.M33;
            m.M34 = a.M34 + b.M34;
            m.M41 = a.M41 + b.M41;
            m.M42 = a.M42 + b.M42;
            m.M43 = a.M43 + b.M43;
            m.M44 = a.M44 + b.M44;
            return m;
        }

        /// <summary>
        /// Adds two matrices together.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Add(ref Matrix a, ref Matrix b, out Matrix result) {
            result.M11 = a.M11 + b.M11;
            result.M12 = a.M12 + b.M12;
            result.M13 = a.M13 + b.M13;
            result.M14 = a.M14 + b.M14;
            result.M21 = a.M21 + b.M21;
            result.M22 = a.M22 + b.M22;
            result.M23 = a.M23 + b.M23;
            result.M24 = a.M24 + b.M24;
            result.M31 = a.M31 + b.M31;
            result.M32 = a.M32 + b.M32;
            result.M33 = a.M33 + b.M33;
            result.M34 = a.M34 + b.M34;
            result.M41 = a.M41 + b.M41;
            result.M42 = a.M42 + b.M42;
            result.M43 = a.M43 + b.M43;
            result.M44 = a.M44 + b.M44;
        }

        /// <summary>
        /// Subtracts one matrix from another.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix Subtract(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = a.M11 - b.M11;
            m.M12 = a.M12 - b.M12;
            m.M13 = a.M13 - b.M13;
            m.M14 = a.M14 - b.M14;
            m.M21 = a.M21 - b.M21;
            m.M22 = a.M22 - b.M22;
            m.M23 = a.M23 - b.M23;
            m.M24 = a.M24 - b.M24;
            m.M31 = a.M31 - b.M31;
            m.M32 = a.M32 - b.M32;
            m.M33 = a.M33 - b.M33;
            m.M34 = a.M34 - b.M34;
            m.M41 = a.M41 - b.M41;
            m.M42 = a.M42 - b.M42;
            m.M43 = a.M43 - b.M43;
            m.M44 = a.M44 - b.M44;
            return m;
        }

        /// <summary>
        /// Subtracts one matrix from another.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Subtract(ref Matrix a, ref Matrix b, out Matrix result) {
            result.M11 = a.M11 - b.M11;
            result.M12 = a.M12 - b.M12;
            result.M13 = a.M13 - b.M13;
            result.M14 = a.M14 - b.M14;
            result.M21 = a.M21 - b.M21;
            result.M22 = a.M22 - b.M22;
            result.M23 = a.M23 - b.M23;
            result.M24 = a.M24 - b.M24;
            result.M31 = a.M31 - b.M31;
            result.M32 = a.M32 - b.M32;
            result.M33 = a.M33 - b.M33;
            result.M34 = a.M34 - b.M34;
            result.M41 = a.M41 - b.M41;
            result.M42 = a.M42 - b.M42;
            result.M43 = a.M43 - b.M43;
            result.M44 = a.M44 - b.M44;
        }

        /// <summary>
        /// Multiplies two matrices together.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>Multiplied matrix</returns>
        public static Matrix Multiply(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = (((a.M11 * b.M11) + (a.M12 * b.M21)) + (a.M13 * b.M31)) + (a.M14 * b.M41);
            m.M12 = (((a.M11 * b.M12) + (a.M12 * b.M22)) + (a.M13 * b.M32)) + (a.M14 * b.M42);
            m.M13 = (((a.M11 * b.M13) + (a.M12 * b.M23)) + (a.M13 * b.M33)) + (a.M14 * b.M43);
            m.M14 = (((a.M11 * b.M14) + (a.M12 * b.M24)) + (a.M13 * b.M34)) + (a.M14 * b.M44);
            m.M21 = (((a.M21 * b.M11) + (a.M22 * b.M21)) + (a.M23 * b.M31)) + (a.M24 * b.M41);
            m.M22 = (((a.M21 * b.M12) + (a.M22 * b.M22)) + (a.M23 * b.M32)) + (a.M24 * b.M42);
            m.M23 = (((a.M21 * b.M13) + (a.M22 * b.M23)) + (a.M23 * b.M33)) + (a.M24 * b.M43);
            m.M24 = (((a.M21 * b.M14) + (a.M22 * b.M24)) + (a.M23 * b.M34)) + (a.M24 * b.M44);
            m.M31 = (((a.M31 * b.M11) + (a.M32 * b.M21)) + (a.M33 * b.M31)) + (a.M34 * b.M41);
            m.M32 = (((a.M31 * b.M12) + (a.M32 * b.M22)) + (a.M33 * b.M32)) + (a.M34 * b.M42);
            m.M33 = (((a.M31 * b.M13) + (a.M32 * b.M23)) + (a.M33 * b.M33)) + (a.M34 * b.M43);
            m.M34 = (((a.M31 * b.M14) + (a.M32 * b.M24)) + (a.M33 * b.M34)) + (a.M34 * b.M44);
            m.M41 = (((a.M41 * b.M11) + (a.M42 * b.M21)) + (a.M43 * b.M31)) + (a.M44 * b.M41);
            m.M42 = (((a.M41 * b.M12) + (a.M42 * b.M22)) + (a.M43 * b.M32)) + (a.M44 * b.M42);
            m.M43 = (((a.M41 * b.M13) + (a.M42 * b.M23)) + (a.M43 * b.M33)) + (a.M44 * b.M43);
            m.M44 = (((a.M41 * b.M14) + (a.M42 * b.M24)) + (a.M43 * b.M34)) + (a.M44 * b.M44);
            return m;
        }

        /// <summary>
        /// Multiplies two matrices together.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Multiply(ref Matrix a, ref Matrix b, out Matrix result) {
            result.M11 = (((a.M11 * b.M11) + (a.M12 * b.M21)) + (a.M13 * b.M31)) + (a.M14 * b.M41);
            result.M12 = (((a.M11 * b.M12) + (a.M12 * b.M22)) + (a.M13 * b.M32)) + (a.M14 * b.M42);
            result.M13 = (((a.M11 * b.M13) + (a.M12 * b.M23)) + (a.M13 * b.M33)) + (a.M14 * b.M43);
            result.M14 = (((a.M11 * b.M14) + (a.M12 * b.M24)) + (a.M13 * b.M34)) + (a.M14 * b.M44);
            result.M21 = (((a.M21 * b.M11) + (a.M22 * b.M21)) + (a.M23 * b.M31)) + (a.M24 * b.M41);
            result.M22 = (((a.M21 * b.M12) + (a.M22 * b.M22)) + (a.M23 * b.M32)) + (a.M24 * b.M42);
            result.M23 = (((a.M21 * b.M13) + (a.M22 * b.M23)) + (a.M23 * b.M33)) + (a.M24 * b.M43);
            result.M24 = (((a.M21 * b.M14) + (a.M22 * b.M24)) + (a.M23 * b.M34)) + (a.M24 * b.M44);
            result.M31 = (((a.M31 * b.M11) + (a.M32 * b.M21)) + (a.M33 * b.M31)) + (a.M34 * b.M41);
            result.M32 = (((a.M31 * b.M12) + (a.M32 * b.M22)) + (a.M33 * b.M32)) + (a.M34 * b.M42);
            result.M33 = (((a.M31 * b.M13) + (a.M32 * b.M23)) + (a.M33 * b.M33)) + (a.M34 * b.M43);
            result.M34 = (((a.M31 * b.M14) + (a.M32 * b.M24)) + (a.M33 * b.M34)) + (a.M34 * b.M44);
            result.M41 = (((a.M41 * b.M11) + (a.M42 * b.M21)) + (a.M43 * b.M31)) + (a.M44 * b.M41);
            result.M42 = (((a.M41 * b.M12) + (a.M42 * b.M22)) + (a.M43 * b.M32)) + (a.M44 * b.M42);
            result.M43 = (((a.M41 * b.M13) + (a.M42 * b.M23)) + (a.M43 * b.M33)) + (a.M44 * b.M43);
            result.M44 = (((a.M41 * b.M14) + (a.M42 * b.M24)) + (a.M43 * b.M34)) + (a.M44 * b.M44);
        }

        /// <summary>
        /// Multiplies a matrix by a scalar value.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="scale">Scalar</param>
        /// <returns>Multiplied matrix</returns>
        public static Matrix Multiply(Matrix value, float scale) {
            Matrix m;
            m.M11 = value.M11 * scale;
            m.M12 = value.M12 * scale;
            m.M13 = value.M13 * scale;
            m.M14 = value.M14 * scale;
            m.M21 = value.M21 * scale;
            m.M22 = value.M22 * scale;
            m.M23 = value.M23 * scale;
            m.M24 = value.M24 * scale;
            m.M31 = value.M31 * scale;
            m.M32 = value.M32 * scale;
            m.M33 = value.M33 * scale;
            m.M34 = value.M34 * scale;
            m.M41 = value.M41 * scale;
            m.M42 = value.M42 * scale;
            m.M43 = value.M43 * scale;
            m.M44 = value.M44 * scale;
            return m;
        }

        /// <summary>
        /// Multiplies a matrix by a scalar value.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="scale">Scalar</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Multiply(ref Matrix value, float scale, out Matrix result) {
            result.M11 = value.M11 * scale;
            result.M12 = value.M12 * scale;
            result.M13 = value.M13 * scale;
            result.M14 = value.M14 * scale;
            result.M21 = value.M21 * scale;
            result.M22 = value.M22 * scale;
            result.M23 = value.M23 * scale;
            result.M24 = value.M24 * scale;
            result.M31 = value.M31 * scale;
            result.M32 = value.M32 * scale;
            result.M33 = value.M33 * scale;
            result.M34 = value.M34 * scale;
            result.M41 = value.M41 * scale;
            result.M42 = value.M42 * scale;
            result.M43 = value.M43 * scale;
            result.M44 = value.M44 * scale;
        }

        /// <summary>
        /// Divides the a matrice's components by the components of another.
        /// </summary>
        /// <param name="a">Source matrix</param>
        /// <param name="b">Divisor</param>
        /// <returns>Divided matrix</returns>
        public static Matrix Divide(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = a.M11 / b.M11;
            m.M12 = a.M12 / b.M12;
            m.M13 = a.M13 / b.M13;
            m.M14 = a.M14 / b.M14;
            m.M21 = a.M21 / b.M21;
            m.M22 = a.M22 / b.M22;
            m.M23 = a.M23 / b.M23;
            m.M24 = a.M24 / b.M24;
            m.M31 = a.M31 / b.M31;
            m.M32 = a.M32 / b.M32;
            m.M33 = a.M33 / b.M33;
            m.M34 = a.M34 / b.M34;
            m.M41 = a.M41 / b.M41;
            m.M42 = a.M42 / b.M42;
            m.M43 = a.M43 / b.M43;
            m.M44 = a.M44 / b.M44;
            return m;
        }

        /// <summary>
        /// Divides the a matrice's components by the components of another.
        /// </summary>
        /// <param name="a">Source matrix</param>
        /// <param name="b">Divisor</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Divide(ref Matrix a, ref Matrix b, out Matrix result) {
            result.M11 = a.M11 / b.M11;
            result.M12 = a.M12 / b.M12;
            result.M13 = a.M13 / b.M13;
            result.M14 = a.M14 / b.M14;
            result.M21 = a.M21 / b.M21;
            result.M22 = a.M22 / b.M22;
            result.M23 = a.M23 / b.M23;
            result.M24 = a.M24 / b.M24;
            result.M31 = a.M31 / b.M31;
            result.M32 = a.M32 / b.M32;
            result.M33 = a.M33 / b.M33;
            result.M34 = a.M34 / b.M34;
            result.M41 = a.M41 / b.M41;
            result.M42 = a.M42 / b.M42;
            result.M43 = a.M43 / b.M43;
            result.M44 = a.M44 / b.M44;
        }

        /// <summary>
        /// Divides the a matrice's components by a scalar value.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="divisor">Divisor</param>
        /// <returns>Divided matrix</returns>
        public static Matrix Divide(Matrix value, float divisor) {
            Matrix m;
            float invScale = 1.0f / divisor;
            m.M11 = value.M11 * invScale;
            m.M12 = value.M12 * invScale;
            m.M13 = value.M13 * invScale;
            m.M14 = value.M14 * invScale;
            m.M21 = value.M21 * invScale;
            m.M22 = value.M22 * invScale;
            m.M23 = value.M23 * invScale;
            m.M24 = value.M24 * invScale;
            m.M31 = value.M31 * invScale;
            m.M32 = value.M32 * invScale;
            m.M33 = value.M33 * invScale;
            m.M34 = value.M34 * invScale;
            m.M41 = value.M41 * invScale;
            m.M42 = value.M42 * invScale;
            m.M43 = value.M43 * invScale;
            m.M44 = value.M44 * invScale;
            return m;
        }

        /// <summary>
        /// Divides the a matrice's components by a scalar value.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="divisor">Divisor</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Divide(ref Matrix value, float divisor, out Matrix result) {
            float invScale = 1.0f / divisor;
            result.M11 = value.M11 * invScale;
            result.M12 = value.M12 * invScale;
            result.M13 = value.M13 * invScale;
            result.M14 = value.M14 * invScale;
            result.M21 = value.M21 * invScale;
            result.M22 = value.M22 * invScale;
            result.M23 = value.M23 * invScale;
            result.M24 = value.M24 * invScale;
            result.M31 = value.M31 * invScale;
            result.M32 = value.M32 * invScale;
            result.M33 = value.M33 * invScale;
            result.M34 = value.M34 * invScale;
            result.M41 = value.M41 * invScale;
            result.M42 = value.M42 * invScale;
            result.M43 = value.M43 * invScale;
            result.M44 = value.M44 * invScale;
        }

        /// <summary>
        /// Decomposes a 3D scale/rotation/translation (SRT) matrix into its
        /// scaling, rotation, and translation components.
        /// </summary>
        /// <param name="scale">Existing vector to hold scale</param>
        /// <param name="rotation">Existing quaternion to hold rotation</param>
        /// <param name="translation">Existing vector to hold translation</param>
        /// <returns>Returns true if rotation was derived, if false the identity quaternion is returned</returns>
        public bool Decompose(out Vector3 scale, out Quaternion rotation, out Vector3 translation) {
            float m11 = this.M11;
            float m12 = this.M12;
            float m13 = this.M13;
            float m14 = this.M14;
            float m21 = this.M21;
            float m22 = this.M22;
            float m23 = this.M23;
            float m24 = this.M24;
            float m31 = this.M31;
            float m32 = this.M32;
            float m33 = this.M33;
            float m34 = this.M34;
            float m41 = this.M41;
            float m42 = this.M42;
            float m43 = this.M43;
            float m44 = this.M44;

            //Flag for if the operation was a success - only returns false if we can't get a rotation
            bool flag = true;

            //Determine which axis to start with dealing with negative scaling
            int index1;
            int index2;
            int index3;

            //Setup a canonical basis
            Basis canonicalBasis = new Basis(Vector3.UnitX, Vector3.UnitY, Vector3.UnitZ);

            //Setup a basis for the rotation matrix
            Vector3 basisX;
            basisX.X = m11;
            basisX.Y = m12;
            basisX.Z = m13;

            Vector3 basisY;
            basisY.X = m21;
            basisY.Y = m22;
            basisY.Z = m23;

            Vector3 basisZ;
            basisZ.X = m31;
            basisZ.Y = m32;
            basisZ.Z = m33;

            Basis rotBasis = new Basis(basisX, basisY, basisZ);

            //Get the translation component
            translation.X = m41;
            translation.Y = m42;
            translation.Z = m43;

            //Get scaling component
            float scaleX = basisX.Length();
            float scaleY = basisY.Length();
            float scaleZ = basisZ.Length();

            if(scaleX < scaleY) {
                if(scaleY < scaleZ) {
                    index1 = 2;
                    index2 = 1;
                    index3 = 0;
                } else {
                    index1 = 1;
                    if(scaleX < scaleZ) {
                        index2 = 2;
                        index3 = 0;
                    } else {
                        index2 = 0;
                        index3 = 2;
                    }
                }
            } else if(scaleX < scaleZ) {
                index1 = 2;
                index2 = 0;
                index3 = 1;
            } else {
                index1 = 0;
                if(scaleY < scaleZ) {
                    index2 = 2;
                    index3 = 1;
                } else {
                    index2 = 1;
                    index3 = 2;
                }
            }
            float temp = 0;
            switch(index1) {
                case 0:
                    temp = scaleX;
                    break;
                case 1:
                    temp = scaleY;
                    break;
                case 2:
                    temp = scaleZ;
                    break;
            }
            bool cb = false;
            if(temp < MathHelper.ZeroTolerance) {
                cb = true;
                rotBasis[index1] = canonicalBasis[index1];
            }
            rotBasis[index1] = Vector3.Normalize(rotBasis[index1]);

            temp = 0;
            switch(index2) {
                case 0:
                    temp = scaleX;
                    break;
                case 1:
                    temp = scaleY;
                    break;
                case 2:
                    temp = scaleZ;
                    break;
            }
            if(temp < MathHelper.ZeroTolerance) {
                int index4;
                float absX = (float) System.Math.Abs((double) rotBasis[index1].X);
                float absY = (float) System.Math.Abs((double) rotBasis[index1].Y);
                float absZ = (float) System.Math.Abs((double) rotBasis[index1].Z);

                if(absX < absY) {
                    if(absY < absZ) {
                        index4 = 0;
                    } else if(absX < absZ) {
                        index4 = 0;
                    } else {
                        index4 = 2;
                    }
                } else if(absX < absZ) {
                    index4 = 1;
                } else if(absY < absZ) {
                    index4 = 1;
                } else {
                    index4 = 2;
                }
                if(cb) {
                    rotBasis[index4] = Vector3.Cross(rotBasis[index2], rotBasis[index1]);
                } 
            }
            rotBasis[index2] = Vector3.Normalize(rotBasis[index2]);

            temp = 0;
            switch(index3) {
                case 0:
                    temp = scaleX;
                    break;
                case 1:
                    temp = scaleY;
                    break;
                case 2:
                    temp = scaleZ;
                    break;
            }
            if(temp < MathHelper.ZeroTolerance) {
                rotBasis[index2] = Vector3.Cross(rotBasis[index3], rotBasis[index1]);
            }
            rotBasis[index3] = Vector3.Normalize(rotBasis[index3]);

            float test = rotBasis.Determinant();
            float det = rotBasis.Determinant();
            if(det < 0.0f) {
                switch(index1) {
                    case 0:
                        scaleX = -scaleX;
                        break;
                    case 1:
                        scaleY = -scaleY;
                        break;
                    case 2:
                        scaleZ = -scaleZ;
                        break;
                }
                rotBasis[index1] = -rotBasis[index1];
                det = -det;
            }
            det--;
            det *= det;
            //Grab rotation
            if(MathHelper.ZeroTolerance < det) {
                rotation = Quaternion.Identity;
                flag = false;
            } else {
                Matrix m = Matrix.Identity;
                m.M11 = rotBasis.XAxis.X;
                m.M12 = rotBasis.XAxis.Y;
                m.M13 = rotBasis.XAxis.Z;
                m.M21 = rotBasis.YAxis.X;
                m.M22 = rotBasis.YAxis.Y;
                m.M23 = rotBasis.YAxis.Z;
                m.M31 = rotBasis.ZAxis.X;
                m.M32 = rotBasis.ZAxis.Y;
                m.M33 = rotBasis.ZAxis.Z;
                Quaternion.FromRotationMatrix(ref m, out rotation);
            }
            //Grab scaling
            scale.X = scaleX;
            scale.Y = scaleY;
            scale.Z = scaleZ;

            return flag;
        }

        /// <summary>
        /// Compute the determinant of this matrix.
        /// </summary>
        /// <returns>Determinant</returns>
        public float Determinant() {
            float m11 = this.M11;
            float m12 = this.M12;
            float m13 = this.M13;
            float m14 = this.M14;
            float m21 = this.M21;
            float m22 = this.M22;
            float m23 = this.M23;
            float m24 = this.M24;
            float m31 = this.M31;
            float m32 = this.M32;
            float m33 = this.M33;
            float m34 = this.M34;
            float m41 = this.M41;
            float m42 = this.M42;
            float m43 = this.M43;
            float m44 = this.M44;

            float h1 = (m33 * m44) - (m34 * m43);
            float h2 = (m32 * m44) - (m34 * m42);
            float h3 = (m32 * m43) - (m33 * m42);
            float h4 = (m31 * m44) - (m34 * m41);
            float h5 = (m31 * m43) - (m33 * m41);
            float h6 = (m31 * m42) - (m32 * m41);
            return ((((m11 * (((m22 * h1) - (m23 * h2)) + (m24 * h3))) - (m12 * (((m21 * h1) - (m23 * h4)) + (m24 * h5)))) + (m13 * (((m21 * h2) - (m22 * h4)) + (m24 * h6)))) - (m14 * (((m21 * h3) - (m22 * h5)) + (m23 * h6))));
        }

        /// <summary>
        /// Compute the inverse of the specified matrix.
        /// </summary>
        /// <param name="matrix">Source matrix</param>
        /// <returns>Inverse of matrix</returns>
        public static Matrix Invert(Matrix matrix) {
            float m11 = matrix.M11;
            float m12 = matrix.M12;
            float m13 = matrix.M13;
            float m14 = matrix.M14;
            float m21 = matrix.M21;
            float m22 = matrix.M22;
            float m23 = matrix.M23;
            float m24 = matrix.M24;
            float m31 = matrix.M31;
            float m32 = matrix.M32;
            float m33 = matrix.M33;
            float m34 = matrix.M34;
            float m41 = matrix.M41;
            float m42 = matrix.M42;
            float m43 = matrix.M43;
            float m44 = matrix.M44;

            float h1 = (m33 * m44) - (m34 * m43);
            float h2 = (m32 * m44) - (m34 * m42);
            float h3 = (m32 * m43) - (m33 * m42);
            float h4 = (m31 * m44) - (m34 * m41);
            float h5 = (m31 * m43) - (m33 * m41);
            float h6 = (m31 * m42) - (m32 * m41);

            float e1 = ((m22 * h1) - (m23 * h2)) + (m24 * h3);
            float e2 = -(((m21 * h1) - (m23 * h4)) + (m24 * h5));
            float e3 = ((m21 * h2) - (m22 * h4)) + (m24 * h6);
            float e4 = -(((m21 * h3) - (m22 * h5)) + (m23 * h6));
            float invDet = 1.0f / ((((m11 * e1) + (m12 * e2)) + (m13 * e3)) + (m14 * e4));

            Matrix m;

            m.M11 = e1 * invDet;
            m.M21 = e2 * invDet;
            m.M31 = e3 * invDet;
            m.M41 = e4 * invDet;

            m.M12 = -(((m12 * h1) - (m13 * h2)) + (m14 * h3)) * invDet;
            m.M22 = (((m11 * h1) - (m13 * h4)) + (m14 * h5)) * invDet;
            m.M32 = -(((m11 * h2) - (m12 * h4)) + (m14 * h6)) * invDet;
            m.M42 = (((m11 * h3) - (m12 * h5)) + (m13 * h6)) * invDet;

            float h7 = (m23 * m44) - (m24 * m43);
            float h8 = (m22 * m44) - (m24 * m42);
            float h9 = (m22 * m43) - (m23 * m42);
            float h10 = (m21 * m44) - (m24 * m41);
            float h11 = (m21 * m43) - (m23 * m41);
            float h12 = (m21 * m42) - (m22 * m41);

            m.M13 = (((m12 * h7) - (m13 * h8)) + (m14 * h9)) * invDet;
            m.M23 = -(((m11 * h7) - (m13 * h10)) + (m14 * h11)) * invDet;
            m.M33 = (((m11 * h8) - (m12 * h10)) + (m14 * h12)) * invDet;
            m.M43 = -(((m11 * h9) - (m12 * h11)) + (m13 * h12)) * invDet;

            float h13 = (m23 * m34) - (m24 * m33);
            float h14 = (m22 * m34) - (m24 * m32);
            float h15 = (m22 * m33) - (m23 * m32);
            float h16 = (m21 * m34) - (m24 * m31);
            float h17 = (m21 * m33) - (m23 * m31);
            float h18 = (m21 * m32) - (m22 * m31);

            m.M14 = -(((m12 * h13) - (m13 * h14)) + (m14 * h15)) * invDet;
            m.M24 = (((m11 * h13) - (m13 * h16)) + (m14 * h17)) * invDet;
            m.M34 = -(((m11 * h14) - (m12 * h16)) + (m14 * h18)) * invDet;
            m.M44 = (((m11 * h15) - (m12 * h17)) + (m13 * h18)) * invDet;
            return m;
        }

        /// <summary>
        /// Compute the inverse of the specified matrix.
        /// </summary>
        /// <param name="matrix">Source matrix</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Invert(ref Matrix matrix, out Matrix result) {
            float m11 = matrix.M11;
            float m12 = matrix.M12;
            float m13 = matrix.M13;
            float m14 = matrix.M14;
            float m21 = matrix.M21;
            float m22 = matrix.M22;
            float m23 = matrix.M23;
            float m24 = matrix.M24;
            float m31 = matrix.M31;
            float m32 = matrix.M32;
            float m33 = matrix.M33;
            float m34 = matrix.M34;
            float m41 = matrix.M41;
            float m42 = matrix.M42;
            float m43 = matrix.M43;
            float m44 = matrix.M44;

            float h1 = (m33 * m44) - (m34 * m43);
            float h2 = (m32 * m44) - (m34 * m42);
            float h3 = (m32 * m43) - (m33 * m42);
            float h4 = (m31 * m44) - (m34 * m41);
            float h5 = (m31 * m43) - (m33 * m41);
            float h6 = (m31 * m42) - (m32 * m41);

            float e1 = ((m22 * h1) - (m23 * h2)) + (m24 * h3);
            float e2 = -(((m21 * h1) - (m23 * h4)) + (m24 * h5));
            float e3 = ((m21 * h2) - (m22 * h4)) + (m24 * h6);
            float e4 = -(((m21 * h3) - (m22 * h5)) + (m23 * h6));
            float invDet = 1.0f / ((((m11 * e1) + (m12 * e2)) + (m13 * e3)) + (m14 * e4));

            result.M11 = e1 * invDet;
            result.M21 = e2 * invDet;
            result.M31 = e3 * invDet;
            result.M41 = e4 * invDet;

            result.M12 = -(((m12 * h1) - (m13 * h2)) + (m14 * h3)) * invDet;
            result.M22 = (((m11 * h1) - (m13 * h4)) + (m14 * h5)) * invDet;
            result.M32 = -(((m11 * h2) - (m12 * h4)) + (m14 * h6)) * invDet;
            result.M42 = (((m11 * h3) - (m12 * h5)) + (m13 * h6)) * invDet;

            float h7 = (m23 * m44) - (m24 * m43);
            float h8 = (m22 * m44) - (m24 * m42);
            float h9 = (m22 * m43) - (m23 * m42);
            float h10 = (m21 * m44) - (m24 * m41);
            float h11 = (m21 * m43) - (m23 * m41);
            float h12 = (m21 * m42) - (m22 * m41);

            result.M13 = (((m12 * h7) - (m13 * h8)) + (m14 * h9)) * invDet;
            result.M23 = -(((m11 * h7) - (m13 * h10)) + (m14 * h11)) * invDet;
            result.M33 = (((m11 * h8) - (m12 * h10)) + (m14 * h12)) * invDet;
            result.M43 = -(((m11 * h9) - (m12 * h11)) + (m13 * h12)) * invDet;

            float h13 = (m23 * m34) - (m24 * m33);
            float h14 = (m22 * m34) - (m24 * m32);
            float h15 = (m22 * m33) - (m23 * m32);
            float h16 = (m21 * m34) - (m24 * m31);
            float h17 = (m21 * m33) - (m23 * m31);
            float h18 = (m21 * m32) - (m22 * m31);

            result.M14 = -(((m12 * h13) - (m13 * h14)) + (m14 * h15)) * invDet;
            result.M24 = (((m11 * h13) - (m13 * h16)) + (m14 * h17)) * invDet;
            result.M34 = -(((m11 * h14) - (m12 * h16)) + (m14 * h18)) * invDet;
            result.M44 = (((m11 * h15) - (m12 * h17)) + (m13 * h18)) * invDet;
        }

        /// <summary>
        /// Linearly interpolate between the two specified matrices.
        /// </summary>
        /// <param name="a">Starting matrix</param>
        /// <param name="b">Ending matrix</param>
        /// <param name="percent">Percent to interpolate by</param>
        /// <returns>Interpolated matrix</returns>
        public static Matrix Lerp(Matrix a, Matrix b, float percent) {
            Matrix m;
            m.M11 = a.M11 + ((b.M11 - a.M11) * percent);
            m.M12 = a.M12 + ((b.M12 - a.M12) * percent);
            m.M13 = a.M13 + ((b.M13 - a.M13) * percent);
            m.M14 = a.M14 + ((b.M14 - a.M14) * percent);
            m.M21 = a.M21 + ((b.M21 - a.M21) * percent);
            m.M22 = a.M22 + ((b.M22 - a.M22) * percent);
            m.M23 = a.M23 + ((b.M23 - a.M23) * percent);
            m.M24 = a.M24 + ((b.M24 - a.M24) * percent);
            m.M31 = a.M31 + ((b.M31 - a.M31) * percent);
            m.M32 = a.M32 + ((b.M32 - a.M32) * percent);
            m.M33 = a.M33 + ((b.M33 - a.M33) * percent);
            m.M34 = a.M34 + ((b.M34 - a.M34) * percent);
            m.M41 = a.M41 + ((b.M41 - a.M41) * percent);
            m.M42 = a.M42 + ((b.M42 - a.M42) * percent);
            m.M43 = a.M43 + ((b.M43 - a.M43) * percent);
            m.M44 = a.M44 + ((b.M44 - a.M44) * percent);
            return m;
        }

        /// <summary>
        /// Linearly interpolate between the two specified matrices.
        /// </summary>
        /// <param name="a">Starting matrix</param>
        /// <param name="b">Ending matrix</param>
        /// <param name="percent">Percent to interpolate by</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Lerp(ref Matrix a, ref Matrix b, float percent, out Matrix result) {
            result.M11 = a.M11 + ((b.M11 - a.M11) * percent);
            result.M12 = a.M12 + ((b.M12 - a.M12) * percent);
            result.M13 = a.M13 + ((b.M13 - a.M13) * percent);
            result.M14 = a.M14 + ((b.M14 - a.M14) * percent);
            result.M21 = a.M21 + ((b.M21 - a.M21) * percent);
            result.M22 = a.M22 + ((b.M22 - a.M22) * percent);
            result.M23 = a.M23 + ((b.M23 - a.M23) * percent);
            result.M24 = a.M24 + ((b.M24 - a.M24) * percent);
            result.M31 = a.M31 + ((b.M31 - a.M31) * percent);
            result.M32 = a.M32 + ((b.M32 - a.M32) * percent);
            result.M33 = a.M33 + ((b.M33 - a.M33) * percent);
            result.M34 = a.M34 + ((b.M34 - a.M34) * percent);
            result.M41 = a.M41 + ((b.M41 - a.M41) * percent);
            result.M42 = a.M42 + ((b.M42 - a.M42) * percent);
            result.M43 = a.M43 + ((b.M43 - a.M43) * percent);
            result.M44 = a.M44 + ((b.M44 - a.M44) * percent);
        }

        /// <summary>
        /// Creates a view matrix with specified eye position, a target, and the up vector in the world.
        /// </summary>
        /// <param name="eyePosition">Position of the camera</param>
        /// <param name="target">Camera's target</param>
        /// <param name="up">Up vector of camera</param>
        /// <returns>View matrix</returns>
        public static Matrix CreateViewMatrix(Vector3 eyePosition, Vector3 target, Vector3 up) {
            Matrix m;
            Vector3 direction = Vector3.Normalize(eyePosition - target);
            Vector3 left = Vector3.Normalize(Vector3.Cross(up, direction));
            Vector3 newUp = Vector3.Cross(direction, left);
            newUp.Normalize();
            m.M11 = left.X;
            m.M12 = newUp.X;
            m.M13 = direction.X;
            m.M14 = 0.0f;
            m.M21 = left.Y;
            m.M22 = newUp.Y;
            m.M23 = direction.Y;
            m.M24 = 0.0f;
            m.M31 = left.Z;
            m.M32 = newUp.Z;
            m.M33 = direction.Z;
            m.M34 = 0.0f;
            m.M41 = -Vector3.Dot(left, eyePosition);
            m.M42 = -Vector3.Dot(newUp, eyePosition);
            m.M43 = -Vector3.Dot(direction, eyePosition);
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a view matrix with specified eye position, a target, and the up vector in the world.
        /// </summary>
        /// <param name="position">Position of the camera</param>
        /// <param name="target">Camera's target</param>
        /// <param name="up">Up vector of camera</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void CreateViewMatrix(ref Vector3 position, ref Vector3 target, ref Vector3 up, out Matrix result) {
            Vector3 direction = Vector3.Normalize(position - target);
            Vector3 left = Vector3.Normalize(Vector3.Cross(up, direction));
            Vector3 newUp = Vector3.Cross(direction, left);
            newUp.Normalize();
            result.M11 = left.X;
            result.M12 = newUp.X;
            result.M13 = direction.X;
            result.M14 = 0.0f;
            result.M21 = left.Y;
            result.M22 = newUp.Y;
            result.M23 = direction.Y;
            result.M24 = 0.0f;
            result.M31 = left.Z;
            result.M32 = newUp.Z;
            result.M33 = direction.Z;
            result.M34 = 0.0f;
            result.M41 = -Vector3.Dot(left, position);
            result.M42 = -Vector3.Dot(newUp, position);
            result.M43 = -Vector3.Dot(direction, position);
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates an orthographic projection matrix.
        /// </summary>
        /// <param name="width">Width of view volume</param>
        /// <param name="height">Height of view volume</param>
        /// <param name="zNearPlane">Minimum z value of the view volume</param>
        /// <param name="zFarPlane">Maximum z value of the view volume</param>
        /// <returns>Resulting projection matrix</returns>
        public static Matrix CreateOrthoMatrix(float width, float height, float zNearPlane, float zFarPlane) {
            Matrix m;
            m.M11 = 2.0f / width;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = 2.0f / height;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = 1.0f / (zNearPlane - zFarPlane);
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = zNearPlane / (zNearPlane - zFarPlane);
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates an orthographic projection matrix.
        /// </summary>
        /// <param name="width">Width of view volume</param>
        /// <param name="height">Height of view volume</param>
        /// <param name="zNearPlane">Minimum z value of the view volume</param>
        /// <param name="zFarPlane">Maximum z value of the view volume</param>
        /// <param name="result">Existing matrix to hold the result</param>
        public static void CreateOrthoMatrix(float width, float height, float zNearPlane, float zFarPlane, out Matrix result) {
            result.M11 = 2.0f / width;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = 2.0f / height;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = 1.0f / (zNearPlane - zFarPlane);
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = zNearPlane / (zNearPlane - zFarPlane);
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a custom orthographic projection matrix.
        /// </summary>
        /// <param name="left">Minimum x value of view volume</param>
        /// <param name="right">Maximum x value of view volume</param>
        /// <param name="bottom">Minimum y value of view volume</param>
        /// <param name="top">Maximum y value of view volume</param>
        /// <param name="zNearPlane">Minimum z value of view volume</param>
        /// <param name="zFarPlane">Maximum z value of view volume</param>
        /// <returns>Resulting projection matrix</returns>
        public static Matrix CreateOrthoMatrix(float left, float right, float bottom, float top, float zNearPlane, float zFarPlane) {
            Matrix m;
            m.M11 = 2.0f / (right - left);
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = 2.0f / (top - bottom);
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = 1.0f / (zNearPlane - zFarPlane);
            m.M34 = 0.0f;
            m.M41 = (left + right) / (left - right);
            m.M42 = (top + bottom) / (bottom - top);
            m.M43 = zNearPlane / (zNearPlane - zFarPlane);
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a custom orthographic projection matrix.
        /// </summary>
        /// <param name="left">Minimum x value of view volume</param>
        /// <param name="right">Maximum x value of view volume</param>
        /// <param name="bottom">Minimum y value of view volume</param>
        /// <param name="top">Maximum y value of view volume</param>
        /// <param name="zNearPlane">Minimum z value of view volume</param>
        /// <param name="zFarPlane">Maximum z value of view volume</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void CreateOrthoMatrix(float left, float right, float bottom, float top, float zNearPlane, float zFarPlane, out Matrix result) {
            result.M11 = 2.0f / (right - left);
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = 2.0f / (top - bottom);
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = 1.0f / (zNearPlane - zFarPlane);
            result.M34 = 0.0f;
            result.M41 = (left + right) / (left - right);
            result.M42 = (top + bottom) / (bottom - top);
            result.M43 = zNearPlane / (zNearPlane - zFarPlane);
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a perspective projection matrix.
        /// </summary>
        /// <param name="width">Width of the view volume at the near view plane</param>
        /// <param name="height">Height of the view volume at the near view plane</param>
        /// <param name="nearPlaneDistance">Distance to the near view plane</param>
        /// <param name="farPlaneDistance">Distance to the far view plane</param>
        /// <returns>Resulting projection matrix</returns>
        public static Matrix CreatePerspectiveMatrix(float width, float height, float nearPlaneDistance, float farPlaneDistance) {
            Matrix m;
            if(nearPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot be negative!");
            }
            if(farPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("farPlaneDistance", "Far plane distance cannot be negative!");
            }
            if(nearPlaneDistance >= farPlaneDistance) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot larger than far plane distance!");
            }
            m.M11 = (2.0f * nearPlaneDistance) / width;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = (2.0f * nearPlaneDistance) / height;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = farPlaneDistance / (nearPlaneDistance - farPlaneDistance);
            m.M34 = -1.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = (nearPlaneDistance * farPlaneDistance) / (nearPlaneDistance - farPlaneDistance);
            m.M44 = 0.0f;
            return m;
        }

        /// <summary>
        /// Creates a perspective projection matrix.
        /// </summary>
        /// <param name="width">Width of the view volume at the near view plane</param>
        /// <param name="height">Height of the view volume at the near view plane</param>
        /// <param name="nearPlaneDistance">Distance to the near view plane</param>
        /// <param name="farPlaneDistance">Distance to the far view plane</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void CreatePerspectiveMatrix(float width, float height, float nearPlaneDistance, float farPlaneDistance, out Matrix result) {
            if(nearPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot be negative!");
            }
            if(farPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("farPlaneDistance", "Far plane distance cannot be negative!");
            }
            if(nearPlaneDistance >= farPlaneDistance) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot larger than far plane distance!");
            }
            result.M11 = (2.0f * nearPlaneDistance) / width;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = (2.0f * nearPlaneDistance) / height;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = farPlaneDistance / (nearPlaneDistance - farPlaneDistance);
            result.M34 = -1.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = (nearPlaneDistance * farPlaneDistance) / (nearPlaneDistance - farPlaneDistance);
            result.M44 = 0.0f;
        }

        /// <summary>
        /// Creates a perspective projection matrix based on a field of view.
        /// </summary>
        /// <param name="fieldOfView">Field of view angle in radians, must be between 0 and Pi</param>
        /// <param name="aspectRatio">Aspect ratio, defined as the view space's width divided by height</param>
        /// <param name="nearPlaneDistance">Distance to the near plane</param>
        /// <param name="farPlaneDistance">Distance to the far plane</param>
        /// <returns>Resulting projection matrix</returns>
        public static Matrix CreatePerspectiveFOVMatrix(float fieldOfView, float aspectRatio, float nearPlaneDistance, float farPlaneDistance) {
            Matrix m;
            if((fieldOfView <= 0.0f) || (fieldOfView >= MathHelper.Pi)) {
                throw new ArgumentOutOfRangeException("fieldOfView", "Field of view out of range, must be between 0 and Pi!");
            }
            if(nearPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot be negative!");
            }
            if(farPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("farPlaneDistance", "Far plane distance cannot be negative!");
            }
            if(nearPlaneDistance >= farPlaneDistance) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot larger than far plane distance!");
            }
            float h = 1.0f / ((float) System.Math.Tan((double) (fieldOfView * 0.5f)));
            float w = h / aspectRatio;

            m.M11 = w;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = h;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = farPlaneDistance / (nearPlaneDistance - farPlaneDistance);
            m.M34 = -1.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = (nearPlaneDistance * farPlaneDistance) / (nearPlaneDistance - farPlaneDistance);
            m.M44 = 0f;
            return m;
        }

        /// <summary>
        /// Creates a perspective projection matrix based on a field of view.
        /// </summary>
        /// <param name="fieldOfView">Field of view angle in radians, must be between 0 and Pi</param>
        /// <param name="aspectRatio">Aspect ratio, defined as the view space's width divided by height</param>
        /// <param name="nearPlaneDistance">Distance to the near plane</param>
        /// <param name="farPlaneDistance">Distance to the far plane</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void CreatePerspectiveFOVMatrix(float fieldOfView, float aspectRatio, float nearPlaneDistance, float farPlaneDistance, out Matrix result) {
            if((fieldOfView <= 0.0f) || (fieldOfView >= MathHelper.Pi)) {
                throw new ArgumentOutOfRangeException("fieldOfView", "Field of view out of range, must be between 0 and Pi!");
            }
            if(nearPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot be negative!");
            }
            if(farPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("farPlaneDistance", "Far plane distance cannot be negative!");
            }
            if(nearPlaneDistance >= farPlaneDistance) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot larger than far plane distance!");
            }
            float h = 1.0f / ((float) System.Math.Tan((double) (fieldOfView * 0.5f)));
            float w = h / aspectRatio;

            result.M11 = w;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = h;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = farPlaneDistance / (nearPlaneDistance - farPlaneDistance);
            result.M34 = -1.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = (nearPlaneDistance * farPlaneDistance) / (nearPlaneDistance - farPlaneDistance);
            result.M44 = 0f;
        }

        /// <summary>
        /// Creates a custom perspective projection matrix.
        /// </summary>
        /// <param name="left">Minimum x value of the view volume at the near view plane</param>
        /// <param name="right">Maximum x value of the view volume at the near view plane</param>
        /// <param name="bottom">Minimum y value of the view volume at the near view plane</param>
        /// <param name="top">Maximum y value of the view volume at the near view plane</param>
        /// <param name="nearPlaneDistance">Distance to the near view plane</param>
        /// <param name="farPlaneDistance">Distance to the far view plane</param>
        /// <returns>Resulting projection matrix</returns>
        public static Matrix CreatePerspectiveMatrix(float left, float right, float bottom, float top, float nearPlaneDistance, float farPlaneDistance) {
            Matrix m;
            if(nearPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot be negative!");
            }
            if(farPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("farPlaneDistance", "Far plane distance cannot be negative!");
            }
            if(nearPlaneDistance >= farPlaneDistance) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot larger than far plane distance!");
            }
            m.M11 = (2.0f * nearPlaneDistance) / (right - left);
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = (2.0f * nearPlaneDistance) / (top - bottom);
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = (left + right) / (right - left);
            m.M32 = (top + bottom) / (top - bottom);
            m.M33 = farPlaneDistance / (nearPlaneDistance - farPlaneDistance);
            m.M34 = -1.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = (nearPlaneDistance * farPlaneDistance) / (nearPlaneDistance - farPlaneDistance);
            m.M44 = 0.0f;
            return m;
        }

        /// <summary>
        /// Creates a custom perspective projection matrix.
        /// </summary>
        /// <param name="left">Minimum x value of the view volume at the near view plane</param>
        /// <param name="right">Maximum x value of the view volume at the near view plane</param>
        /// <param name="bottom">Minimum y value of the view volume at the near view plane</param>
        /// <param name="top">Maximum y value of the view volume at the near view plane</param>
        /// <param name="nearPlaneDistance">Distance to the near view plane</param>
        /// <param name="farPlaneDistance">Distance to the far view plane</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void CreatePerspectiveMatrix(float left, float right, float bottom, float top, float nearPlaneDistance, float farPlaneDistance, out Matrix result) {
            if(nearPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot be negative!");
            }
            if(farPlaneDistance <= 0.0f) {
                throw new ArgumentOutOfRangeException("farPlaneDistance", "Far plane distance cannot be negative!");
            }
            if(nearPlaneDistance >= farPlaneDistance) {
                throw new ArgumentOutOfRangeException("nearPlaneDistance", "Near plane distance cannot larger than far plane distance!");
            }
            result.M11 = (2.0f * nearPlaneDistance) / (right - left);
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = (2.0f * nearPlaneDistance) / (top - bottom);
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = (left + right) / (right - left);
            result.M32 = (top + bottom) / (top - bottom);
            result.M33 = farPlaneDistance / (nearPlaneDistance - farPlaneDistance);
            result.M34 = -1.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = (nearPlaneDistance * farPlaneDistance) / (nearPlaneDistance - farPlaneDistance);
            result.M44 = 0.0f;
        }

        /// <summary>
        /// Creates a rotation matrix where the object is facing the target along its z axis.
        /// If your object's "forward" facing is down -Z, then the object will correctly face the target.
        /// </summary>
        /// <param name="position">Position of object</param>
        /// <param name="target">Position of target</param>
        /// <param name="worldUp">World yAxis vector</param>
        /// <returns>Resulting rotation matrix</returns>
        public static Matrix LookAt(Vector3 position, Vector3 target, Vector3 worldUp) {
            Vector3 zAxis;
            Vector3 xAxis;
            Vector3 yAxis;

            Vector3.Subtract(ref position, ref target, out zAxis);
            zAxis.Normalize();

            Vector3 worldUpXdir;
            Vector3.Cross(ref worldUp, ref zAxis, out worldUpXdir);
            Vector3.Normalize(ref worldUpXdir, out xAxis);
            Vector3.Cross(ref zAxis, ref xAxis, out yAxis);

            Matrix m;
            FromAxes(ref xAxis, ref yAxis, ref zAxis, out m);
            return m;
        }

        /// <summary>
        /// Creates a rotation matrix where the object is facing the target along its z axis.
        /// If your object's "forward" facing is down -Z, then the object will correctly face the target.
        /// </summary>
        /// <param name="position">Position of object</param>
        /// <param name="target">Position of target</param>
        /// <param name="worldUp">World yAxis vector</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void LookAt(ref Vector3 position, ref Vector3 target, ref Vector3 worldUp, out Matrix result){
            Vector3 zAxis;
            Vector3 xAxis;
            Vector3 yAxis;

            Vector3.Subtract(ref position, ref target, out zAxis);
            zAxis.Normalize();

            Vector3 worldUpXdir;
            Vector3.Cross(ref worldUp, ref zAxis, out worldUpXdir);
            Vector3.Normalize(ref worldUpXdir, out xAxis);
            Vector3.Cross(ref zAxis, ref xAxis, out yAxis);

            FromAxes(ref xAxis, ref yAxis, ref zAxis, out result);
        }

        /// <summary>
        /// Creates a rotation matrix from 3 orthogonal axes.
        /// </summary>
        /// <param name="xAxis">X axis (first row, right)</param>
        /// <param name="yAxis">Y axis (second row, up)</param>
        /// <param name="zAxis">Z axis (third row, backward)</param>
        /// <returns>Rotation matrix</returns>
        public static Matrix FromAxes(Vector3 xAxis, Vector3 yAxis, Vector3 zAxis) {
            Matrix m;
            //Set x axis to xAxis (first row)
            m.M11 = xAxis.X;
            m.M12 = xAxis.Y;
            m.M13 = xAxis.Z;
            m.M14 = 0.0f;

            //Set y axis to yAxis (second row)
            m.M21 = yAxis.X;
            m.M22 = yAxis.Y;
            m.M23 = yAxis.Z;
            m.M24 = 0.0f;

            //Set z axis to zAxis (third row)
            m.M31 = zAxis.X;
            m.M32 = zAxis.Y;
            m.M33 = zAxis.Z;
            m.M34 = 0.0f;

            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;

            return m;
        }

        /// <summary>
        /// Creates a rotation matrix from 3 orthogonal axes.
        /// </summary>
        /// <param name="xAxis">X axis (first row, right)</param>
        /// <param name="yAxis">Y axis (second row, up)</param>
        /// <param name="zAxis">Z axis (third row, backward)</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromAxes(ref Vector3 xAxis, ref Vector3 yAxis, ref Vector3 zAxis, out Matrix result) {
            //Set x axis to xAxis (first row)
            result.M11 = xAxis.X;
            result.M12 = xAxis.Y;
            result.M13 = xAxis.Z;
            result.M14 = 0.0f;

            //Set y axis to yAxis (second row)
            result.M21 = yAxis.X;
            result.M22 = yAxis.Y;
            result.M23 = yAxis.Z;
            result.M24 = 0.0f;

            //Set z axis to zAxis (third row)
            result.M31 = zAxis.X;
            result.M32 = zAxis.Y;
            result.M33 = zAxis.Z;
            result.M34 = 0.0f;

            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Transforms the supplied matrix by a quaternion rotation.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="rot">Quaternion rotation</param>
        /// <returns>Transformed matrix</returns>
        public static Matrix Transform(Matrix value, Quaternion rot) {
            float x2 = rot.X + rot.X;
            float y2 = rot.Y + rot.Y;
            float z2 = rot.Z + rot.Z;

            float wx2 = rot.W * x2;
            float wy2 = rot.W * y2;
            float wz2 = rot.W * z2;

            float xx2 = rot.X * x2;
            float xy2 = rot.X * y2;
            float xz2 = rot.X * z2;

            float yy2 = rot.Y * y2;
            float yz2 = rot.Y * z2;

            float zz2 = rot.Z * z2;

            float h1 = (1.0f - yy2) - zz2;
            float h2 = xy2 - wz2;
            float h3 = xz2 + wy2;
            float h4 = xy2 + wz2;
            float h5 = (1.0f - xx2) - zz2;
            float h6 = yz2 - wx2;
            float h7 = xz2 - wy2;
            float h8 = yz2 + wx2;
            float h9 = (1.0f - xx2) - yy2;

            Matrix m;
            m.M11 = ((value.M11 * h1) + (value.M12 * h2)) + (value.M13 * h3);
            m.M12 = ((value.M11 * h4) + (value.M12 * h5)) + (value.M13 * h6);
            m.M13 = ((value.M11 * h7) + (value.M12 * h8)) + (value.M13 * h9);
            m.M14 = value.M14;
            m.M21 = ((value.M21 * h1) + (value.M22 * h2)) + (value.M23 * h3);
            m.M22 = ((value.M21 * h4) + (value.M22 * h5)) + (value.M23 * h6);
            m.M23 = ((value.M21 * h7) + (value.M22 * h8)) + (value.M23 * h9);
            m.M24 = value.M24;
            m.M31 = ((value.M31 * h1) + (value.M32 * h2)) + (value.M33 * h3);
            m.M32 = ((value.M31 * h4) + (value.M32 * h5)) + (value.M33 * h6);
            m.M33 = ((value.M31 * h7) + (value.M32 * h8)) + (value.M33 * h9);
            m.M34 = value.M34;
            m.M41 = ((value.M41 * h1) + (value.M42 * h2)) + (value.M43 * h3);
            m.M42 = ((value.M41 * h4) + (value.M42 * h5)) + (value.M43 * h6);
            m.M43 = ((value.M41 * h7) + (value.M42 * h8)) + (value.M43 * h9);
            m.M44 = value.M44;
            return m;
        }

        /// <summary>
        /// Transforms the supplied matrix by a quaternion rotation.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="rot">Quaternion rotation</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Transform(ref Matrix value, ref Quaternion rot, out Matrix result) {
            float x2 = rot.X + rot.X;
            float y2 = rot.Y + rot.Y;
            float z2 = rot.Z + rot.Z;

            float wx2 = rot.W * x2;
            float wy2 = rot.W * y2;
            float wz2 = rot.W * z2;

            float xx2 = rot.X * x2;
            float xy2 = rot.X * y2;
            float xz2 = rot.X * z2;

            float yy2 = rot.Y * y2;
            float yz2 = rot.Y * z2;

            float zz2 = rot.Z * z2;

            float h1 = (1f - yy2) - zz2;
            float h2 = xy2 - wz2;
            float h3 = xz2 + wy2;
            float h4 = xy2 + wz2;
            float h5 = (1f - xx2) - zz2;
            float h6 = yz2 - wx2;
            float h7 = xz2 - wy2;
            float h8 = yz2 + wx2;
            float h9 = (1f - xx2) - yy2;

            result.M11 = ((value.M11 * h1) + (value.M12 * h2)) + (value.M13 * h3);
            result.M12 = ((value.M11 * h4) + (value.M12 * h5)) + (value.M13 * h6);
            result.M13 = ((value.M11 * h7) + (value.M12 * h8)) + (value.M13 * h9);
            result.M14 = value.M14;
            result.M21 = ((value.M21 * h1) + (value.M22 * h2)) + (value.M23 * h3);
            result.M22 = ((value.M21 * h4) + (value.M22 * h5)) + (value.M23 * h6);
            result.M23 = ((value.M21 * h7) + (value.M22 * h8)) + (value.M23 * h9);
            result.M24 = value.M24;
            result.M31 = ((value.M31 * h1) + (value.M32 * h2)) + (value.M33 * h3);
            result.M32 = ((value.M31 * h4) + (value.M32 * h5)) + (value.M33 * h6);
            result.M33 = ((value.M31 * h7) + (value.M32 * h8)) + (value.M33 * h9);
            result.M34 = value.M34;
            result.M41 = ((value.M41 * h1) + (value.M42 * h2)) + (value.M43 * h3);
            result.M42 = ((value.M41 * h4) + (value.M42 * h5)) + (value.M43 * h6);
            result.M43 = ((value.M41 * h7) + (value.M42 * h8)) + (value.M43 * h9);
            result.M44 = value.M44;
        }

        /// <summary>
        /// Swaps the rows and columns of the supplied matrix.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <returns>Transposed matrix</returns>
        public static Matrix Transpose(Matrix value) {
            Matrix m;
            m.M11 = value.M11;
            m.M12 = value.M21;
            m.M13 = value.M31;
            m.M14 = value.M41;
            m.M21 = value.M12;
            m.M22 = value.M22;
            m.M23 = value.M32;
            m.M24 = value.M42;
            m.M31 = value.M13;
            m.M32 = value.M23;
            m.M33 = value.M33;
            m.M34 = value.M43;
            m.M41 = value.M14;
            m.M42 = value.M24;
            m.M43 = value.M34;
            m.M44 = value.M44;
            return m;
        }

        /// <summary>
        /// Swaps the rows and columns of the supplied matrix.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Transpose(ref Matrix value, out Matrix result) {
            result.M11 = value.M11;
            result.M12 = value.M21;
            result.M13 = value.M31;
            result.M14 = value.M41;
            result.M21 = value.M12;
            result.M22 = value.M22;
            result.M23 = value.M32;
            result.M24 = value.M42;
            result.M31 = value.M13;
            result.M32 = value.M23;
            result.M33 = value.M33;
            result.M34 = value.M43;
            result.M41 = value.M14;
            result.M42 = value.M24;
            result.M43 = value.M34;
            result.M44 = value.M44;
        }

        /// <summary>
        /// Flips the signs of each matrix element.
        /// </summary>
        public void Negate() {
            M11 = -M11;
            M12 = -M12;
            M13 = -M13;
            M14 = -M14;
            M21 = -M21;
            M22 = -M22;
            M23 = -M23;
            M24 = -M24;
            M31 = -M31;
            M32 = -M32;
            M33 = -M33;
            M34 = -M34;
            M41 = -M41;
            M42 = -M42;
            M43 = -M43;
            M44 = -M44;
        }

        /// <summary>
        /// Flips the signs of each element of the supplied
        /// matrix.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <returns>Negated matrix</returns>
        public static Matrix Negate(Matrix value) {
            Matrix m;
            m.M11 = -value.M11;
            m.M12 = -value.M12;
            m.M13 = -value.M13;
            m.M14 = -value.M14;
            m.M21 = -value.M21;
            m.M22 = -value.M22;
            m.M23 = -value.M23;
            m.M24 = -value.M24;
            m.M31 = -value.M31;
            m.M32 = -value.M32;
            m.M33 = -value.M33;
            m.M34 = -value.M34;
            m.M41 = -value.M41;
            m.M42 = -value.M42;
            m.M43 = -value.M43;
            m.M44 = -value.M44;
            return m;
        }

        /// <summary>
        /// Flips the signs of each element of the supplied
        /// matrix.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void Negate(ref Matrix value, out Matrix result) {
            result.M11 = -value.M11;
            result.M12 = -value.M12;
            result.M13 = -value.M13;
            result.M14 = -value.M14;
            result.M21 = -value.M21;
            result.M22 = -value.M22;
            result.M23 = -value.M23;
            result.M24 = -value.M24;
            result.M31 = -value.M31;
            result.M32 = -value.M32;
            result.M33 = -value.M33;
            result.M34 = -value.M34;
            result.M41 = -value.M41;
            result.M42 = -value.M42;
            result.M43 = -value.M43;
            result.M44 = -value.M44;
        }

        /// <summary>
        /// Creates a matrix that represents a rotation around the x-axis.
        /// </summary>
        /// <param name="radians">Angle to rotate, in radians</param>
        /// <returns>Resulting rotation matrix</returns>
        public static Matrix FromRotationX(float radians) {
            Matrix m;
            float cos = (float) System.Math.Cos((double) radians);
            float sin = (float) System.Math.Sin((double) radians);
            m.M11 = 1.0f;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = cos;
            m.M23 = sin;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = -sin;
            m.M33 = cos;
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a matrix that represents a rotation around the x-axis.
        /// </summary>
        /// <param name="radians">Angle to rotate, in radians</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromRotationX(float radians, out Matrix result) {
            float cos = (float) System.Math.Cos((double) radians);
            float sin = (float) System.Math.Sin((double) radians);
            result.M11 = 1.0f;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = cos;
            result.M23 = sin;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = -sin;
            result.M33 = cos;
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a matrix that represents a rotation around the y-axis.
        /// </summary>
        /// <param name="radians">Angle to rotate, in radians</param>
        /// <returns>Resulting rotation matrix</returns>
        public static Matrix FromRotationY(float radians) {
            Matrix m;
            float cos = (float) System.Math.Cos((double) radians);
            float sin = (float) System.Math.Sin((double) radians);
            m.M11 = cos;
            m.M12 = 0.0f;
            m.M13 = -sin;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = 1.0f;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = sin;
            m.M32 = 0.0f;
            m.M33 = cos;
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a matrix that represents a rotation around the y-axis.
        /// </summary>
        /// <param name="radians">Angle to rotate, in radians</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromRotationY(float radians, out Matrix result) {
            float cos = (float) System.Math.Cos((double) radians);
            float sin = (float) System.Math.Sin((double) radians);
            result.M11 = cos;
            result.M12 = 0.0f;
            result.M13 = -sin;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = 1.0f;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = sin;
            result.M32 = 0.0f;
            result.M33 = cos;
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a matrix that represents a rotation around the z-axis.
        /// </summary>
        /// <param name="radians">Angle to rotate, in radians</param>
        /// <returns>Resulting rotation matrix</returns>
        public static Matrix FromRotationZ(float radians) {
            Matrix m;
            float cos = (float) System.Math.Cos((double) radians);
            float sin = (float) System.Math.Sin((double) radians);
            m.M11 = cos;
            m.M12 = sin;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = -sin;
            m.M22 = cos;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = 1.0f;
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a matrix that represents a rotation around the z-axis.
        /// </summary>
        /// <param name="radians">Angle to rotate, in radians</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromRotationZ(float radians, out Matrix result) {
            float cos = (float) System.Math.Cos((double) radians);
            float sin = (float) System.Math.Sin((double) radians);
            result.M11 = cos;
            result.M12 = sin;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = -sin;
            result.M22 = cos;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = 1.0f;
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a scaling matrix from the supplied scale vector.
        /// </summary>
        /// <param name="scale">Scaling vector where each component corresponds to each axis to scale on</param>
        /// <returns>Resulting scaling matrix</returns>
        public static Matrix FromScale(Vector3 scale) {
            Matrix m;
            float x = scale.X;
            float y = scale.Y;
            float z = scale.Z;
            m.M11 = x;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = y;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = z;
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a scaling matrix from the supplied scale vector.
        /// </summary>
        /// <param name="scale">Scaling vector where each component corresponds to each axis to scale on</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromScale(ref Vector3 scale, out Matrix result) {
            float x = scale.X;
            float y = scale.Y;
            float z = scale.Z;
            result.M11 = x;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = y;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = z;
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a scaling matrix from the supplied scaling value for uniform scale.
        /// </summary>
        /// <param name="scale">Scaling value for x,y,z axes</param>
        /// <returns>Scaling matrix</returns>
        public static Matrix FromScale(float scale) {
            Matrix m;
            m.M11 = scale;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = scale;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = scale;
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a scaling matrix from the supplied scaling value for uniform scale.
        /// </summary>
        /// <param name="scale">Scaling value for x,y,z axes</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromScale(float scale, out Matrix result) {
            result.M11 = scale;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = scale;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = scale;
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a scale matrix from the supplied scalars.
        /// </summary>
        /// <param name="x">Amount to scale along x axis</param>
        /// <param name="y">Amount to scale along y axis</param>
        /// <param name="z">Amount to scale along z axis</param>
        /// <returns>Resulting scale matrix</returns>
        public static Matrix FromScale(float x, float y, float z) {
            Matrix m;
            m.M11 = x;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = y;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = z;
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a scale matrix from the supplied scalars.
        /// </summary>
        /// <param name="x">Amount to scale along x axis</param>
        /// <param name="y">Amount to scale along y axis</param>
        /// <param name="z">Amount to scale along z axis</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromScale(float x, float y, float z, out Matrix result) {
            result.M11 = x;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = y;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = z;
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a translation matrix from the supplied vector.
        /// </summary>
        /// <param name="translation">Translation vector</param>
        /// <returns>Resulting translation matrix</returns>
        public static Matrix FromTranslation(Vector3 translation) {
            Matrix m;
            m.M11 = 1.0f;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = 1.0f;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = 1.0f;
            m.M34 = 0.0f;
            m.M41 = translation.X;
            m.M42 = translation.Y;
            m.M43 = translation.Z;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a translation matrix from the supplied vector.
        /// </summary>
        /// <param name="translation">Translation vector</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromTranslation(ref Vector3 translation, out Matrix result) {
            result.M11 = 1.0f;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = 1.0f;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = 1.0f;
            result.M34 = 0.0f;
            result.M41 = translation.X;
            result.M42 = translation.Y;
            result.M43 = translation.Z;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a translation matrix from the supplied values.
        /// </summary>
        /// <param name="x">Amount to translate along x axis</param>
        /// <param name="y">Amount to translate along y axis</param>
        /// <param name="z">Amount to translate along z axis</param>
        /// <returns>Resulting translation matrix</returns>
        public static Matrix FromTranslation(float x, float y, float z) {
            Matrix m;
            m.M11 = 1.0f;
            m.M12 = 0.0f;
            m.M13 = 0.0f;
            m.M14 = 0.0f;
            m.M21 = 0.0f;
            m.M22 = 1.0f;
            m.M23 = 0.0f;
            m.M24 = 0.0f;
            m.M31 = 0.0f;
            m.M32 = 0.0f;
            m.M33 = 1.0f;
            m.M34 = 0.0f;
            m.M41 = x;
            m.M42 = y;
            m.M43 = z;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a translation matrix from the supplied values.
        /// </summary>
        /// <param name="x">Amount to translate along x axis</param>
        /// <param name="y">Amount to translate along y axis</param>
        /// <param name="z">Amount to translate along z axis</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromTranslation(float x, float y, float z, out Matrix result) {
            result.M11 = 1.0f;
            result.M12 = 0.0f;
            result.M13 = 0.0f;
            result.M14 = 0.0f;
            result.M21 = 0.0f;
            result.M22 = 1.0f;
            result.M23 = 0.0f;
            result.M24 = 0.0f;
            result.M31 = 0.0f;
            result.M32 = 0.0f;
            result.M33 = 1.0f;
            result.M34 = 0.0f;
            result.M41 = x;
            result.M42 = y;
            result.M43 = z;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a rotation matrix from the supplied quaternion.
        /// </summary>
        /// <param name="rot">Rotation quaternion</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix FromQuaternion(Quaternion rot) {
            float x = rot.X;
            float y = rot.Y;
            float z = rot.Z;
            float w = rot.W;

            float xx = x * x;
            float xy = x * y;
            float xw = x * w;

            float yy = y * y;
            float yw = y * w;
            float yz = y * z;

            float zx = z * x;
            float zz = z * z;
            float zw = z * w;

            Matrix m;
            m.M11 = 1.0f - (2.0f * (yy + zz));
            m.M12 = 2.0f * (xy + zw);
            m.M13 = 2.0f * (zx - yw);
            m.M14 = 0.0f;

            m.M21 = 2.0f * (xy - zw);
            m.M22 = 1.0f - (2.0f * (zz + xx));
            m.M23 = 2.0f * (yz + xw);
            m.M24 = 0.0f;

            m.M31 = 2.0f * (zx + yw);
            m.M32 = 2.0f * (yz - xw);
            m.M33 = 1.0f - (2.0f * (yy + xx));
            m.M34 = 0.0f;

            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a rotation matrix from the supplied quaternion.
        /// </summary>
        /// <param name="rot">Rotation quaternion</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromQuaternion(ref Quaternion rot, out Matrix result) {
            float norm = (rot.W * rot.W) + (rot.X * rot.X) + (rot.Y * rot.Y) + (rot.Z * rot.Z);

            float s = 0;

            //We must ensure the quaternion is normalized first. We also use 2 / norm to save mults later
            if(norm == 1.0f) {
                s = 2;
            } else if(norm > 0.0f) {
                s = 2 / norm;
            }

            //Compute xs/ys/zs since we use them 2-4 times, saves 6 mults
            float xs = rot.X * s;
            float ys = rot.Y * s;
            float zs = rot.Z * s;

            float xx = rot.X * xs;
            float xy = rot.X * ys;
            float xz = rot.X * zs;
            float xw = rot.W * xs;

            float yy = rot.Y * ys;
            float yz = rot.Y * zs;
            float yw = rot.W * ys;

            float zz = rot.Z * zs;
            float zw = rot.W * zs;

            result.M11 = 1.0f - (yy + zz);
            result.M12 =        (xy + zw);
            result.M13 =        (xz - yw);
            result.M14 = 0.0f;

            result.M21 =        (xy - zw);
            result.M22 = 1.0f - (zz + xx);
            result.M23 =        (yz + xw);
            result.M24 = 0.0f;

            result.M31 =        (xz + yw);
            result.M32 =        (yz - xw);
            result.M33 = 1.0f - (yy + xx);
            result.M34 = 0.0f;

            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a rotation matrix from a supplied angle in radians and
        /// axis to rotate about.
        /// </summary>
        /// <param name="angle">Angle in radians</param>
        /// <param name="axis">Axis to rotate about</param>
        /// <returns>Rotation matrix</returns>
        public static Matrix FromAngleAxis(float angle, Vector3 axis) {
            Matrix m;
            float x = axis.X;
            float y = axis.Y;
            float z = axis.Z;

            float sin = (float) System.Math.Sin((double) angle);
            float cos = (float) System.Math.Cos((double) angle);

            float xx = x * x;
            float yy = y * y;
            float zz = z * z;
            float xy = x * y;
            float xz = x * z;
            float yz = y * z;

            m.M11 = xx + (cos * (1.0f - xx));
            m.M12 = (xy - (cos * xy)) + (sin * z);
            m.M13 = (xz - (cos * xz)) - (sin * y);
            m.M14 = 0.0f;
            m.M21 = (xy - (cos * xy)) - (sin * z);
            m.M22 = yy + (cos * (1.0f - yy));
            m.M23 = (yz - (cos * yz)) + (sin * x);
            m.M24 = 0.0f;
            m.M31 = (xz - (cos * xz)) + (sin * y);
            m.M32 = (yz - (cos * yz)) - (sin * x);
            m.M33 = zz + (cos * (1.0f - zz));
            m.M34 = 0.0f;
            m.M41 = 0.0f;
            m.M42 = 0.0f;
            m.M43 = 0.0f;
            m.M44 = 1.0f;
            return m;
        }

        /// <summary>
        /// Creates a rotation matrix from a supplied angle in radians and
        /// axis to rotate about.
        /// </summary>
        /// <param name="angle">Angle in radians</param>
        /// <param name="axis">Axis to rotate about</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromAngleAxis(float angle, ref Vector3 axis, out Matrix result) {
            float x = axis.X;
            float y = axis.Y;
            float z = axis.Z;

            float sin = (float) System.Math.Sin((double) angle);
            float cos = (float) System.Math.Cos((double) angle);

            float xx = x * x;
            float yy = y * y;
            float zz = z * z;
            float xy = x * y;
            float xz = x * z;
            float yz = y * z;

            result.M11 = xx + (cos * (1.0f - xx));
            result.M12 = (xy - (cos * xy)) + (sin * z);
            result.M13 = (xz - (cos * xz)) - (sin * y);
            result.M14 = 0.0f;
            result.M21 = (xy - (cos * xy)) - (sin * z);
            result.M22 = yy + (cos * (1.0f - yy));
            result.M23 = (yz - (cos * yz)) + (sin * x);
            result.M24 = 0.0f;
            result.M31 = (xz - (cos * xz)) + (sin * y);
            result.M32 = (yz - (cos * yz)) - (sin * x);
            result.M33 = zz + (cos * (1.0f - zz));
            result.M34 = 0.0f;
            result.M41 = 0.0f;
            result.M42 = 0.0f;
            result.M43 = 0.0f;
            result.M44 = 1.0f;
        }

        /// <summary>
        /// Creates a rotation matrix from the corresponding yaw, pitch, and roll
        /// angles. The order of the angles are applied in that order.
        /// </summary>
        /// <param name="yaw">Angle to rotate about the y-axis</param>
        /// <param name="pitch">Angle to rotate about the x-axis</param>
        /// <param name="roll">Angle to rotate about the z-axis</param>
        /// <returns>Rotation matrix</returns>
        public static Matrix FromEulerAngles(float yaw, float pitch, float roll) {
            Matrix matrix;
            Quaternion quaternion;
            Quaternion.FromEulerAngles(yaw, pitch, roll, out quaternion);
            FromQuaternion(ref quaternion, out matrix);
            return matrix;
        }

        /// <summary>
        /// Creates a rotation matrix from the corresponding yaw, pitch, and roll
        /// angles. The order of the angles are applied in that order.
        /// </summary>
        /// <param name="yaw">Angle to rotate about the y-axis</param>
        /// <param name="pitch">Angle to rotate about the x-axis</param>
        /// <param name="roll">Angle to rotate about the z-axis</param>
        /// <param name="result">Existing matrix to hold result</param>
        public static void FromEulerAngles(float yaw, float pitch, float roll, out Matrix result) {
            Quaternion quaternion;
            Quaternion.FromEulerAngles(yaw, pitch, roll, out quaternion);
            FromQuaternion(ref quaternion, out result);
        }

        /// <summary>
        /// Tests equality with another matrix.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <returns>True if equal</returns>
        public bool Equals(Matrix value) {
            return (((this.M11 == value.M11) && (this.M12 == value.M12) && (this.M13 == value.M13) && (this.M14 == value.M14))
                && ((this.M21 == value.M21) && (this.M22 == value.M22) && (this.M23 == value.M23) && (this.M24 == value.M24))
                && ((this.M31 == value.M31) && (this.M32 == value.M32) && (this.M33 == value.M33) && (this.M34 == value.M34))
                && ((this.M41 == value.M41) && (this.M42 == value.M42) && (this.M43 == value.M43) && (this.M44 == value.M44)));
        }

        /// <summary>
        /// Tests equality between two matrices.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>True if equal</returns>
        public static bool operator==(Matrix a, Matrix b){
            return (((a.M11 == b.M11) && (a.M12 == b.M12) && (a.M13 == b.M13) && (a.M14 == b.M14))
                && ((a.M21 == b.M21) && (a.M22 == b.M22) && (a.M23 == b.M23) && (a.M24 == b.M24))
                && ((a.M31 == b.M31) && (a.M32 == b.M32) && (a.M33 == b.M33) && (a.M34 == b.M34))
                && ((a.M41 == b.M41) && (a.M42 == b.M42) && (a.M43 == b.M43) && (a.M44 == b.M44)));
        }

        /// <summary>
        /// Tests inequality between two matrices.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>True if not equal</returns>
        public static bool operator!=(Matrix a, Matrix b) {
            return !(a == b);
        }

        /// <summary>
        /// Adds two matrices together.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator+(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = a.M11 + b.M11;
            m.M12 = a.M12 + b.M12;
            m.M13 = a.M13 + b.M13;
            m.M14 = a.M14 + b.M14;
            m.M21 = a.M21 + b.M21;
            m.M22 = a.M22 + b.M22;
            m.M23 = a.M23 + b.M23;
            m.M24 = a.M24 + b.M24;
            m.M31 = a.M31 + b.M31;
            m.M32 = a.M32 + b.M32;
            m.M33 = a.M33 + b.M33;
            m.M34 = a.M34 + b.M34;
            m.M41 = a.M41 + b.M41;
            m.M42 = a.M42 + b.M42;
            m.M43 = a.M43 + b.M43;
            m.M44 = a.M44 + b.M44;
            return m;
        }

        /// <summary>
        /// Flips the signs of the components of the specified matrix.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator-(Matrix value) {
            Matrix m;
            m.M11 = -value.M11;
            m.M12 = -value.M12;
            m.M13 = -value.M13;
            m.M14 = -value.M14;
            m.M21 = -value.M21;
            m.M22 = -value.M22;
            m.M23 = -value.M23;
            m.M24 = -value.M24;
            m.M31 = -value.M31;
            m.M32 = -value.M32;
            m.M33 = -value.M33;
            m.M34 = -value.M34;
            m.M41 = -value.M41;
            m.M42 = -value.M42;
            m.M43 = -value.M43;
            m.M44 = -value.M44;
            return m;
        }

        /// <summary>
        /// Subtracts a matrix from the other.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator-(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = a.M11 - b.M11;
            m.M12 = a.M12 - b.M12;
            m.M13 = a.M13 - b.M13;
            m.M14 = a.M14 - b.M14;
            m.M21 = a.M21 - b.M21;
            m.M22 = a.M22 - b.M22;
            m.M23 = a.M23 - b.M23;
            m.M24 = a.M24 - b.M24;
            m.M31 = a.M31 - b.M31;
            m.M32 = a.M32 - b.M32;
            m.M33 = a.M33 - b.M33;
            m.M34 = a.M34 - b.M34;
            m.M41 = a.M41 - b.M41;
            m.M42 = a.M42 - b.M42;
            m.M43 = a.M43 - b.M43;
            m.M44 = a.M44 - b.M44;
            return m;
        }

        /// <summary>
        /// Multiplies two matrices.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator*(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = (((a.M11 * b.M11) + (a.M12 * b.M21)) + (a.M13 * b.M31)) + (a.M14 * b.M41);
            m.M12 = (((a.M11 * b.M12) + (a.M12 * b.M22)) + (a.M13 * b.M32)) + (a.M14 * b.M42);
            m.M13 = (((a.M11 * b.M13) + (a.M12 * b.M23)) + (a.M13 * b.M33)) + (a.M14 * b.M43);
            m.M14 = (((a.M11 * b.M14) + (a.M12 * b.M24)) + (a.M13 * b.M34)) + (a.M14 * b.M44);
            m.M21 = (((a.M21 * b.M11) + (a.M22 * b.M21)) + (a.M23 * b.M31)) + (a.M24 * b.M41);
            m.M22 = (((a.M21 * b.M12) + (a.M22 * b.M22)) + (a.M23 * b.M32)) + (a.M24 * b.M42);
            m.M23 = (((a.M21 * b.M13) + (a.M22 * b.M23)) + (a.M23 * b.M33)) + (a.M24 * b.M43);
            m.M24 = (((a.M21 * b.M14) + (a.M22 * b.M24)) + (a.M23 * b.M34)) + (a.M24 * b.M44);
            m.M31 = (((a.M31 * b.M11) + (a.M32 * b.M21)) + (a.M33 * b.M31)) + (a.M34 * b.M41);
            m.M32 = (((a.M31 * b.M12) + (a.M32 * b.M22)) + (a.M33 * b.M32)) + (a.M34 * b.M42);
            m.M33 = (((a.M31 * b.M13) + (a.M32 * b.M23)) + (a.M33 * b.M33)) + (a.M34 * b.M43);
            m.M34 = (((a.M31 * b.M14) + (a.M32 * b.M24)) + (a.M33 * b.M34)) + (a.M34 * b.M44);
            m.M41 = (((a.M41 * b.M11) + (a.M42 * b.M21)) + (a.M43 * b.M31)) + (a.M44 * b.M41);
            m.M42 = (((a.M41 * b.M12) + (a.M42 * b.M22)) + (a.M43 * b.M32)) + (a.M44 * b.M42);
            m.M43 = (((a.M41 * b.M13) + (a.M42 * b.M23)) + (a.M43 * b.M33)) + (a.M44 * b.M43);
            m.M44 = (((a.M41 * b.M14) + (a.M42 * b.M24)) + (a.M43 * b.M34)) + (a.M44 * b.M44);
            return m;
        }

        /// <summary>
        /// Multiplies a matrix by a scalar value.
        /// </summary>
        /// <param name="value">Source matrix</param>
        /// <param name="scale">Scalar</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator*(Matrix value, float scale) {
            Matrix m;
            m.M11 = value.M11 * scale;
            m.M12 = value.M12 * scale;
            m.M13 = value.M13 * scale;
            m.M14 = value.M14 * scale;
            m.M21 = value.M21 * scale;
            m.M22 = value.M22 * scale;
            m.M23 = value.M23 * scale;
            m.M24 = value.M24 * scale;
            m.M31 = value.M31 * scale;
            m.M32 = value.M32 * scale;
            m.M33 = value.M33 * scale;
            m.M34 = value.M34 * scale;
            m.M41 = value.M41 * scale;
            m.M42 = value.M42 * scale;
            m.M43 = value.M43 * scale;
            m.M44 = value.M44 * scale;
            return m;
        }

        /// <summary>
        /// Multiplies a matrix by a scalar value.
        /// </summary>
        /// <param name="scale">Scalar</param>
        /// <param name="value">Source matrix</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator*(float scale, Matrix value) {
            Matrix m;
            m.M11 = value.M11 * scale;
            m.M12 = value.M12 * scale;
            m.M13 = value.M13 * scale;
            m.M14 = value.M14 * scale;
            m.M21 = value.M21 * scale;
            m.M22 = value.M22 * scale;
            m.M23 = value.M23 * scale;
            m.M24 = value.M24 * scale;
            m.M31 = value.M31 * scale;
            m.M32 = value.M32 * scale;
            m.M33 = value.M33 * scale;
            m.M34 = value.M34 * scale;
            m.M41 = value.M41 * scale;
            m.M42 = value.M42 * scale;
            m.M43 = value.M43 * scale;
            m.M44 = value.M44 * scale;
            return m;
        }

        /// <summary>
        /// Divides the components of a matrix by the counterpart components
        /// of the other.
        /// </summary>
        /// <param name="a">First matrix</param>
        /// <param name="b">Second matrix (divisor)</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator/(Matrix a, Matrix b) {
            Matrix m;
            m.M11 = a.M11 / b.M11;
            m.M12 = a.M12 / b.M12;
            m.M13 = a.M13 / b.M13;
            m.M14 = a.M14 / b.M14;
            m.M21 = a.M21 / b.M21;
            m.M22 = a.M22 / b.M22;
            m.M23 = a.M23 / b.M23;
            m.M24 = a.M24 / b.M24;
            m.M31 = a.M31 / b.M31;
            m.M32 = a.M32 / b.M32;
            m.M33 = a.M33 / b.M33;
            m.M34 = a.M34 / b.M34;
            m.M41 = a.M41 / b.M41;
            m.M42 = a.M42 / b.M42;
            m.M43 = a.M43 / b.M43;
            m.M44 = a.M44 / b.M44;
            return m;
        }

        /// <summary>
        /// Divides a matrix by a scalar.
        /// </summary>
        /// <param name="a">Source matrix</param>
        /// <param name="divisor">Divisor</param>
        /// <returns>Resulting matrix</returns>
        public static Matrix operator/(Matrix a, float divisor) {
            float invDivisor = 1.0f / divisor;
            Matrix m;
            m.M11 = a.M11 * invDivisor;
            m.M12 = a.M12 * invDivisor;
            m.M13 = a.M13 * invDivisor;
            m.M14 = a.M14 * invDivisor;
            m.M21 = a.M21 * invDivisor;
            m.M22 = a.M22 * invDivisor;
            m.M23 = a.M23 * invDivisor;
            m.M24 = a.M24 * invDivisor;
            m.M31 = a.M31 * invDivisor;
            m.M32 = a.M32 * invDivisor;
            m.M33 = a.M33 * invDivisor;
            m.M34 = a.M34 * invDivisor;
            m.M41 = a.M41 * invDivisor;
            m.M42 = a.M42 * invDivisor;
            m.M43 = a.M43 * invDivisor;
            m.M44 = a.M44 * invDivisor;
            return m;
        }

        /// <summary>
        /// Tests equality between this matrix and the specified object.
        /// </summary>
        /// <param name="obj">Object to compare</param>
        /// <returns>True if equal</returns>
        public override bool Equals(object obj) {
            if(obj is Matrix) {
                return Equals((Matrix) obj);
            }
            return false;
        }

        /// <summary>
        /// Get the hash code of this object.
        /// </summary>
        /// <returns>Hash code</returns>
        public override int GetHashCode() {
            return (((((((((((((((this.M11.GetHashCode() + this.M12.GetHashCode()) + this.M13.GetHashCode()) + this.M14.GetHashCode())
                + this.M21.GetHashCode()) + this.M22.GetHashCode()) + this.M23.GetHashCode()) + this.M24.GetHashCode())
                + this.M31.GetHashCode()) + this.M32.GetHashCode()) + this.M33.GetHashCode()) + this.M34.GetHashCode())
                + this.M41.GetHashCode()) + this.M42.GetHashCode()) + this.M43.GetHashCode()) + this.M44.GetHashCode());
        }

        /// <summary>
        /// Get the string representation of this object.
        /// </summary>
        /// <returns>String representation</returns>
        public override String ToString() {
            CultureInfo info = CultureInfo.CurrentCulture;
            Object[] args = new object[] { this.M11.ToString(info), this.M12.ToString(info), this.M13.ToString(info), this.M14.ToString(info), this.M21.ToString(info), this.M22.ToString(info), this.M23.ToString(info), this.M24.ToString(info), this.M31.ToString(info), this.M32.ToString(info), this.M33.ToString(info), this.M34.ToString(info), this.M41.ToString(info), this.M42.ToString(info), this.M43.ToString(info), this.M44.ToString(info) };
            return String.Format(info, "{{[M11:{0} M12:{1} M13:{2} M14:{3}] [M21:{4} M22:{5} M23:{6} M24:{7}] [M31:{8} M32:{9} M33:{10} M34:{11}] [M41:{12} M42:{13} M43:{14} M44:{15}]}}", args);
        }
    }

    //Basis struct for use in Matrix.Decompose()
    internal struct Basis {
        public Vector3 XAxis;
        public Vector3 YAxis;
        public Vector3 ZAxis;

        public Vector3 this[int index] {
            get {
                switch(index) {
                    case 0:
                        return XAxis;
                    case 1:
                        return YAxis;
                    case 2:
                        return ZAxis;
                    default:
                        return Vector3.Zero;
                }
            }
            set {
                switch(index) {
                    case 0:
                        XAxis = value;
                        break;
                    case 1:
                        YAxis = value;
                        break;
                    case 2:
                        ZAxis = value;
                        break;
                }
            }
        }

        public Basis(Vector3 v1, Vector3 v2, Vector3 v3) {
            XAxis = v1;
            YAxis = v2;
            ZAxis = v3;
        }

        public float Determinant() {
            float m11 = XAxis.X;
            float m12 = XAxis.Y;
            float m13 = XAxis.Z;
            float m21 = YAxis.X;
            float m22 = YAxis.Y;
            float m23 = YAxis.Z;
            float m31 = ZAxis.X;
            float m32 = ZAxis.Y;
            float m33 = ZAxis.Z;

            float h1 = m33;
            float h2 = m32;
            float h4 = m31;
            return ((((m11 * (((m22 * h1) - (m23 * h2)))) - (m12 * (((m21 * h1) - (m23 * h4))))) + (m13 * (((m21 * h2) - (m22 * h4))))));
        }
    }
}
